NASA Technical Memorandum 113112 
ICOMP-97-08; CMOTT-97-03 


Modeling of Turbulent Swirling Flows 


Tsan-Hsing Shih, Jiang Zhu, and William Liou 
Institute for Computational Mechanics in Propulsion 
and Center for Modeling of Turbulence and Transition 
Cleveland, Ohio 

Kuo-Huey Chen 
University of Toledo 
Toldeo, Ohio 

Nan-Suey Liu 
Lewis Research Center 
Cleveland, Ohio 

John L. Lumley 
Cornell University 
Ithaca, New York 


August 1997 



Space Administration 



Modeling of turbulent swirling flows 


Tsan-Hsing Shih, Jiang Zhu and William Liou 

Center for Modeling of Turbulence and Transition , ICOMP, NASA Lewis Research Center , 22800 

Cedar Point Rd Brook Park , OH 4414% 

Kuo-Huey Chen 

The University of Toledo , Ohio 


Nan-Suey Liu 

NASA Lewis Research Center , Cleveland , OH 4413$ 

John L. Lumley 

Cornell University, Ithaca , New York 


Abstract 

Aircraft engine combustors generally involve turbulent swirling flows in order to enhance fuel-air 
mixing and flame stabilization. It has long been recognized that eddy viscosity turbulence models 
are unable to appropriately model swirling flows. Therefore, it has been suggested that, for the 
modeling of these flows, a second order closure scheme should be considered because of its ability 
in the modeling of rotational and curvature effects. However, this scheme will require solution 
of many complicated second moment transport equations (six Reynolds stresses plus other scalar 
fluxes and variances), which is a difficult task for any CFD implementations. Also, this scheme will 
require a large amount of computer resources for a general combustor swirling flow. 

This report is devoted to the development of a cubic Reynolds stress-strain model for turbulent 
swirling flows, and was inspired by the work of Launder’s group at UMIST. Using this type of 
model, one only needs to solve two turbulence equations, one for the turbulent kinetic energy k 
and the other for the dissipation rate e. The cubic model developed in this report is based on a 
general Reynolds stress-strain relationship (Shih and Lumley, 1993). Two flows have been chosen 
for model evaluation. One is a fully developed rotating pipe flow, and the other is a more complex 
flow with swirl and recirculation. 
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1 Introduction 

For better fuel-air mixing and flame stabilization in a combustor, a swirl is generally asso- 
ciated with the flows. Therefore, accurate modeling of turbulent swirling flows is important 
in engine combustor design. Common turbulence models used in engineering calculations 
are eddy viscosity models which include zero-equation and two-equation models (e.g., mix- 
ing length models and k-e models). However, it has long been recognized that this type of 
eddy viscosity model is not appropriate for predicting swirling flows. In fact, the deficiency 
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of eddy viscosity models for swirling flows can be analytically demonstrated by modeling 
a fully developed rotating pipe flow (Fu, 1995). Measured swirl velocity in the pipe varies 
approximately as the square of the normalized radius (r 2 ), however, eddy viscosity models 
produce an exact linear profile of the swirl velocity, which describes a solid body rotation. 

To avoid this kind of deficiency of eddy viscosity models, a second order closure scheme 
has been suggested for modeling of swirling flows because of its ability to simulate the 
effects of mean rotation and curvature. However, this requires solving many complicated 
second moment transport equations, which involve six Reynolds stresses plus other scalar 
fluxes and variances. Because of this complexity and because of the large computer resources 
required, second moment transport equation models have not been successfully implemented 
in combustor swirling flows. 

Recent developments in nonlinear Reynolds stress-strain models bring a practical method 
for combustion flow calculations because of their potential in simulating turbulent swirling 
flows with only two modeled turbulence transport equations (Craft et al, 1993). Further 
development and evaluation of these models are of great interest to both CFD development 
and modern aircraft engine combustor design. 

The model developed in this report is based on a general Reynolds stress-strain relationship 
which is an explicit expression for the Reynolds stresses in terms of a tensorial polynomial 
of mean velocity gradients. It is derived from a generalized Cayley-Hamilton relation. This 
general formulation contains terms up to the sixth power of the mean velocity gradient with 
eleven undetermined coefficients. Obviously, for any practical application, we need to trun- 
cate this polynomial. Shih, Zhu and Lumley (1995) suggested a quadratic formulation and 
determined the three relevant coefficients by using the realizability constraints of Reynolds 
stresses and a result from rapid distortion theory analysis. This quadratic model works quite 
successfully for many complex flows including flows with separation. However, our recent 
calculations of swirling flows show that the swirl velocity is not appropriately predicted, 
which verifies the finding from Launder’s group at UMIST. Launder (1995) pointed out that 
“the weaknesses of the linear eddy viscosity model can not be rectified by introducing just 
quadratic terms to the stress-strain relation.” 

In this report, we retain the cubic terms from a general Reynolds stress-stain formulation 
and determine the coefficients by using a similar method used in Shih et al’s quadratic model 
and the measured data from rotating pipe flows. Modeled k-e equations are used together 
with the cubic Reynolds stress-strain model for mean flow calculations. The first test flow is 
that of fully developed pipe flow rotating about its own axial axis with various rotation rates 
(Imao, Itoh and Harada, 1996). The second test flow is a more complex flow with swirl and 
recirculation (Roback and Johnson, 1983). These two flows both have detailed experimental 
data on mean velocity components. The comparisons between the experimental data and 
computational results from models will be reported in detail. 

In this report, there axe four appendices. In Appendix A, the derivation of the proposed cubic 
model is described. Appendix B gives the equations in a general coordinate system, which 
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will be useful for studying flows in various curvilinear coordinate systems. For example, 
axisymmetric flows will be most conveniently studied in a cylindrical coordinate system. 
Therefore, in Appendix C and Appendix D, we write the equations for a general flow and 
an axisymmetric flow respectively in a cylindrical coordinate system. 


2 Cubic Reynolds stress model 


A cubic Reynolds stress model, used in this study for modeling of turbulent swirling flows, 
is developed in Appendix A. The resultant cubic model can be expressed in terms of mean 
velocity gradients, [Ty, or in terms of mean strain and rotation rates, and Q,ij. Here, we 
list both forms for convenience of their applications. 

In terms of mean velocity gradients, the cubic model for Reynolds stresses is 
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where II S is defined in Eq. (12). 

In terms of mean strain and rotation rates, Eq. (1) can be written as 
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I Is is the second principal invariant of S^- defined as 
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Note that in the above equations, S k k means Sii + S 2 2 + S 3 3 and Sl k means S lp S pl + S 2p S p2 + 
S 3 P Spz in which each term contains a summation operator on the subscript “p” . 

It should also be mentioned that the eddy viscosity p T in Eq. (3) will become the standard 

k 2 

form of px = pC „ — for high turbulent Reynolds number flows ( — >> 1). 
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3 Modeling of turbulent swirling flows 


The model proposed in the previous section will be used for modeling of swirling flows in 
this study. The first flow is a fully developed rotating pipe flow (Imao, Itoh and Harada, 
1996). This flow was used for model development; however, a pipe flow with various axial 
rotating rates is still a critical test case for the model. The second flow is a more complex 
swirling flow with recirculation and separation (Roback and Johnson, 1983), which is often 
encountered in an aircraft engine combustor. 


3.1 Rotating pipe flow 


A fully developed rotating pipe flow provides a very clean test case for checking the turbu- 
lence model’s ability to model swirling flows. As mentioned previously, commonly used eddy 
viscosity models fail to predict this flow. In fact, one can show that any eddy viscosity model 
will produce a solution of solid body rotation for a rotating pipe flow, while experimental 
data shows that the flow is not a solid body rotation. Experiments further demonstrate 
that the characteristics of a pipe flow changes significantly with the axial rotation rate. For 
example, for a fixed mass flux, the axial rotation will strongly reduce the pressure drop. In 
other words, for a fixed pressure drop, the axial rotation will increase the total mass flux. 
However, standard eddy viscosity models show no such changes at all. 

In a fully developed turbulent pipe flow, all the axial gradients, d/dx , and the azimuthal 
derivatives, d/d9, are zero, and so is the radial velocity V = 0. The non-zero velocity 
components are the axial velocity U and the tangential (or swirl) velocity W = rfi, where 
Q, is the angular velocity. Equations for this flow are 
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where 5 — ^2SijSij = y (^) + (^r — - 7 -) • The nonlinear parts of turbulent stresses, r xr 
and To r , from the proposed cubic model, Eq. (1) or Eq. (10), are 


T XT — 0 (1^) 



The production rate of turbulent kinetic energy P k is 
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= py/kyjp. Other model constants 


used in this report are standard: a k = 1 , a t = 1.3, C el = 1.44 and C t2 = 1.92. Depending on 
particular modeled k-s equations, the model coefficients and damping function f ^ may have 
different formulations proposed by various researchers. For example, if Eq. (17) (Shih et al, 
1995) is used together with 
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From Eq. (14), it is easy to show that any eddy viscosity model will produce a solution of 
solid body rotation, i.e., W/W wa u = r/R, where W wa u is the swirl velocity of the wall and R 
is the radius of the pipe. It can also be shown that any quadratic Reynolds stress models will 
have no contributions to the component t$ t for a fully develed rotating pipe flow. Therefore, 
they will also produce a solution of solid body rotation, just like an eddy viscosity model 
does. Equations (13)-(17) can be easily and accurately solved by a parabolic code. Figures 
1-3 show the results of the present cubic model with Eqs. (15) and (17) compared with 
the measurements by Imao, et al (1996). The results from the standard k-e eddy viscosity 
model are also included for comparison. In the figures, the rotation parameter N is defined 
a s N = W wa u/U m , where U m is the average velocity of the pipe. The Reynolds number based 
on U m and R is 20000. As shown in these figures, the standard k-e eddy viscosity model has 
totally missed the effect of axial rotations on the pipe flow. In contrast, the present cubic 
Reynolds stress model can capture all the effects of the axial rotation on the pipe flow: it 
increases the centerline velocity and changes the axial velocity profile towards a parabolic 
shape, it maintains non-solid body swirl velocity profile, and it reduces the relative turbulent 
kinetic energy kjU\. 


3.2 Complex swirling flow with recirculation 


A confined swirling coaxial jet was experimentally studied by Roback and Johnson (1983). 
Figure 4 shows the general features of the flow. At the inlet, an inner jet and an annular jet 
are ejected into an enlarged duct. Besides an annular separation due to sudden expansion of 
the duct, a central recirculation bubble is created by the swirling flow. This flow feature is 
often observed in an aircraft engine combustor. In this figure, calculated velocity vectors in 
an axisymmetric plane from the cubic model is compared with the one from the standard k-e 
eddy viscosity model. Solutions were obtained by two Navier-Stokes codes. One is CORSAIR 
(Liu et al, 1996) and the other is FAST-2D (Zhu, 1991). Eq. (16) and Eq. (17) are respectively 
used in this calculation. Numerical results from the two codes are quite close to each other. 
Figure 5 compares the calculations of the centerline velocity using a standard k-e eddy 
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viscosity model (SKE) and the present cubic model with the experimental data. The negative 
velocity indicates the central recirculation. It is seen that both models predict the strength 
of central recirculation quite well, but the present model predicts the rear stagnation point 
much better than does the SKE model. This is also reflected in Fig. 4 that the recirculation 
bubble from the cubic model is larger than that from the standard SKE model. Figure 6 shows 
the comparison of calculated and measured mean velocity profiles at x=51mm. Both models 
give reasonably good profiles which are within experimental scatter. However, significant 
differences in the tangential velocity profile between the two models have been found in the 
downstream region. For example, Fig. 7 shows the swirl velocity profile at x=305mm. SKE 
model predicts a nearly solid body rotation, whereas the cubic model shows a non-solid body 
rotation which is consistent with experimental observation. 


4 Conclusion and discussion 


This study shows that nonlinear cubic Reynolds stress-strain models with modeled k-e equa- 
tions have the potential to simulate turbulent swirling flows encountered in aircraft engine 
combustors. The model proposed in this report appears simple and numerically robust in 
CFD applications in which the aircraft engine industry is particularly interested. However, 
further evaluations against other flows are needed in order to determine the flow range of 
the model’s validity and to seek possible further improvements. 

The proposed cubic Reynolds stress model can be combined with existing k-e model equa- 
tions, yet the best combination needs further studies and evaluations. 

The proposed cubic model appears the simplest among other cubic or higher order models; 
however it requires about 15% more CPU time than does a linear k-e eddy viscosity model for 
a general 2D axisymmetric swirling flow. We expect that if a higher order model (e.g., fourth 
or fifth order) is used, then the CPU time for calculating Reynolds stresses will significantly 
increase and the model may become very costly for the calculation of a general 3D swirling 
flow. 
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Fig. 2. Tangential velocity profile in a rotating pipe 
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(b). SKE model 


Fig. 4. Velocity vectors in an axisymmetric plane, (a) from present model, (b) from SKE model. 
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A Appendix: Development of a Cubic Turbulent Model 


A truncated general cubic turbulent stress-strain relation from Shih and Lumley (1993) can 
be written as 


2 k 2 ( 2 \ 

-puiuj = - -pk Sij + CuP— [Uij + Uj yi - -U kyk 6ijJ 
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The six model coefficients in Eq.(A.l) will be determined by the following procedure. First, 
we consider two extreme cases: a pure strain flow and a pure shear flow, and apply realizabil- 
ity constraints on the Reynolds stresses to ensure positive energy components and Schwarz’ 
inequality. This was suggested by Reynolds (1987) and Shih et al (1995), which will allow us 
to determine the model coefficients of C C\, Ci and C 3 . The second procedure is to deter- 
mine the model coefficients C4 and C 5 by using the experimental data of a fully developed 
rotating pipe flow. To analize the pure strain and pure shear flows, it is more convenient to 
write Eq.(A.l) in terms of mean strain and rotation rates, as in the following: 
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(A.3) 
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Note that 5 *■ , SjJ*\ and S- 3 * 5 are all traceless tensors. Using Cayley-Hamilton relation, 
5j - J 5 Sj + 7/sS, - = 0 (A-4) 

S\j*^ can be expressed in terms of quadratic and linear terms as 
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Using Eq.(A.5), we may write Eq.(A.2) as 
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A\ = -(2Ci + ^2 + C3) , A 2 = — (2Ci — C 2 — C3) 

A 3 = C 3 -C 2 , A 4 = (C 5 -C 4 ), A 5 = (C4 + C7 5 ) (A. 8) 

A result from a rapid distortion theory analysis (Reynolds, 1987) states that isotropic tur- 
bulence should not be affected by a pure mean rotation. To satisfy this result, the simplest 
way is to eliminate the pure rotation term in Eq.(A.7), i.e., A 2 = 0, which indicates that 
2Ci = C 2 + C 3 . 

To determine the model coefficients, let us first consider a pure strain flow, in which Ct,j = 0. 
Under this situation, 

fm&J = \pk6ij - (CV - ^IJ s A 5 )^2S*. - (A 1 + -I s A 5 )^-2Slp (A.9) 

In principal axes of S*j , we may write (see Shih, Zhu and Lumley, 1995) 
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where a and b can take on arbitrary values. Then, one may write 
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Following Shill, Zhu and Lumley (1995), for simplicity we set (3 = 0, i.e., A x 

jj* 

which indicates C 2 + C 3 = -~I S A 5 . Then, a must be less than unity, i.e., 

£ 

(CV - — II s A 5 ) < 
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(1995), or Reynolds (1987)]: 
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and can be calculated using the following relations [see Shih et al 
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From Eq.(A.19), can be written as 


C a < 


» - 4 kS 

A s 


k2 rr a 
T H — tIIsAs 


(A.21) 


Now, let us consider a pure shear flow, in which there is only one non-zero component, Uip, 
i.e., 
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In this case, 5 i2 = fl 12 = -U h2 . Under this situation, we obtain from Eq.(A.7) 
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Note that in Eqs. (A.22)-(A.24), the condition (Ai + -/sA 5 ) = 0 has been used, and note also 
that A3 must be positive since the shear Ui f2 will make u\ increase and u\ decrease. Applying 
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Schwarz’ inequality, (U1U2) 2 < u \ u 2 > t° the above equations, we obtain a constraint for A3: 
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Noting that ( 5 *) 2 = 2S'i2-Si2 and Q*S* = 2fi 12 5i2 for the pure shear flow, a generalized 
expression for A 3 may be written as 
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To ensure a positive real value of A 3 , the coefficient C M must be also restricted by the 
following condition for any large values of 5*: 


V2 





(A. 28) 


The formulations for C M and A 3 , i.e., Eqs. (A.21) and (A.26), will ensure realizability of 
turbulent stresses. However, A 4 and A5 are left to be further determined, which are related 
to the coefficients C4 and C5 by Eq.(A.8). 

To determine A 4 and A 5 , or C 4 and C 5 , let us study a fully developed rotating pipe flow. 
In this case, only two components of the non-linear part of turbulent stresses, r xr and T 0r . 
appear in the mean flow equations, i.e., Eqs. (13) and (14), which are 


T*- — C+ 


pk 4 W (dW W\du 


f dW _ w\ 

\ dr r ) 


dr 


(A. 29) 


^ pk* Tzr dW fdW W' 
t 0t = - C A ~-W- 

£ 6 or \ or r 


-C, 


pk 4 


dr y dr 
2 

1 mi \ 

W 


duy xir dw fdw 

dr ) dr 1 dr 




(A. 30) 


Now integrate the Eq. (14) for the velocity W component at a steady state to obtain 
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(A.31) 


(H + h t )t 




-c 5 



aw (aw 
w — — — — 

dr y dr 



W ( r \ 2 

Experimental data show that — — « ( — ) for a large range of W wall . Here, R is the radius 

W wa u \R / 

of the pipe, W wa u is the wall swirl velocity. Insert this relation into the above equation, we 
obtain, for high turbulent Reynolds numbers, 


Rt 


pk 4 


T 2 Ri WLu ~ 


pk 4 


r 2 

2—W, 


R 4 


wall 


+ 




o 


(A. 32) 


If we write 


pk 4 


2jjlWw,=a’M r, C. 


pk 4 


2 jpWl all + 


(& 


— a' ^t 


(A.33) 


then from Eq.(A.32), we must require a 1 + (3' ~ 1. The coefficients C\ and C5 can be expressed 
as 


, _ sL 

e 3 R 4 


P' t*T 


pk 4 

[2^7" + (s 

rf 

£ 3 

R 4 \c 

b J 


(A.34) 


In a fully developed, rotating, pipe flow, we find that the following relations hold, 


9 rjwLu 

R 4 


q* q* 


flij j 


(auV 

R 1 \9r) 


= i (7SJSJ + n^i) 


(A. 35) 
(A. 36) 


Finally, we obtain expressions for C 4 and C5 as follows 


Mr 

frj |(s -) 2 - (si *) 2 

P' 

^ ( 7 <s-) 2 + (n*) 2 ) 


(A.37) 
(A. 38) 


19 



From the calculation of rotating pipe flows, we find that the following coefficients seem 
appropriate (i.e., we set a' = 0, (3' = 1.6): 


C - 1 

M “ kU* 

(A. 39) 

4.0 + A s 

P 


) 

O 

II 

O" 

(A.40) 

r 1-6 Pt 

5 ~ pk 4 7(S*) 2 + (n*) 2 

(A.41) 


e 3 4 

where 

U* = ^5*5* + = fl* = (A.42) 

Equations (A.40) and (A. 8) suggest that A4 = A5 = C5. 

Now, we may summarize the cubic model and its coefficients as follows: 

2 pJc 2 

pU{Uj = -pk8{j + pr2, S+j + A 3 — 

o 

-2A 5 ^(n ik s 2 kj - s? k n kj + n ik s km a mj 

-^klSl^mkSij + IIsS *•) (A.43) 

where 




n 

1 

^ ” , n A kU* 
4.0 + A s 

€ 

Al .\ 

i_Wi£y 

2 M ( £ J 


0.5+^fTS* 
2 e 2 

A — _ 

1.6px 

— 

pk 4 7(S*) 2 + (fl*) 2 


4 


(A.44) 

(A.45) 


(A.46) 

(A.47) 


In Eq. (A.43), we have used the fact that S* k Cl k j — Sli k S£j = Si k Q. k j — Q ik S k j and — 

= QikSlj — Sf k Ct k j. In addition, C M must also be constrained by the conditions from 
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(A. 48) 


Eqs. (A.21) and (A.28), i.e., 


^ < \ As 


kS*' 


-i 


k 2 

+ -^UsAs 


j ^ V2 fkS* 

and C M < — I 

o V S 


) 


-1 


The cubic model can be directly expressed in terms of mean velocity gradients, i.e., Eq. 
(A.l). The corresponding coefficients are 


C i = - 

c 2 = - 


lk r a 

2l IsA > 


- ( ^3 H Is A 5 



then the cubic model, Eq.(A.43), becomes 


2 / 2 \ 
-puiuj = - -pk 6jj + p T [Uij + Ujj - -U k ,k 

Mf,e (UuVkj-Utj.Uj^ 


+ 


2 £ 2 



U k ,iU k , p U P j + U k jU k , p U v , 



-\h {VijVH + Uij.Ut, , - |n,««) 

-i/ s (mmO« + - §iw«)] 


% 


(A. 49) 
(A. 50) 

(A.51) 

(A. 52) 
(A. 53) 


(A.54) 
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B Appendix: Equations in a General Coordinate System 


In this appendix, a set of mean flow equations with a general cubic model will be written in a 
general coordinate system. This appendix will be found useful for studying turbulent flows in 
a curvilinear coordinate system. We start with the governing equations in general tensorial 
form. In Appendix C, we will write these equations in a cylindrical coordinate system as an 
example to show how to write the equations and models for a specific curvilinear coordinate 
system. 


B.l Equations in tensorial form 


e,, + (pU‘).=o 


where 


P,i + 9 ir 

P + Uj,i g P 

(B.l) 

(B.2) 


) kj\ + P k -pe 

(B.3) 

A a k/ 

’ r 

0* + £) 

£ i] +C cl h e -P k -C t2 f 2 p £ - 



+C £3 ^5 jV 5 J 5, r 

(B.4) 

s = J 

'2 StjS'j 

(B-5) 


The turbulent stress is written in the following form: 

-py.u, = - 2 -pkg XJ + fir (Uij + U si - + T ij 


(B.6) 


where the subscript “ t ,” denotes a tensorial derivative, and gij are the two metric tensors 
of a coordinate system, which are defined in Eq.(B.16). The nonlinear part of the general 
cubic model, r^, is 


T<i = + u iM u$ - jn, 9ij ) 

+ C ^ (g u U i j t U i j - in 2 
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+ c 4 ^ + g-'u^uiu,, - In, g 

+ c s ^ (up^u? + vp^u* - ?n 5 Siy ) 

(B.7) 

where, 


n< = 0JD} , n 2 = cf'upn , n 3 = g u up iim ux 

(B.8) 

In addition, the often used scalar parameters S* and fl* defined in 
Eq.(A.20) can be written as 

Eq.(A.42) and W* in 

(S ’) 2 = 5 + Up$ - i((/‘) 2 

(B.9) 

( n -) 2 = i («*>'%, 

(B.10) 

o* o* c* 

TT7. Jj „kl „TTin 

W* = g 3 g g 

(B.ll) 


The nonlinear part of turbulent stress r tJ -, Eq. (B.7), can also be expressed in terms of mean 
strain and rotation rates 5,j and Clij which will be listed in Eq. (B.33). 


B.2 Equations in a general coodinate system 


Let x l represent a general curvilinear coordinate system, then the corresponding contravari- 

dx* 

ant velocity is defined as U x = — and the covariant velocity is defined as £/, = g^U 3 . To 

dt 

write Eqs. (B.l)-(B.IO) in this general coordinate system, we need the following expressions 
for various tensorial derivatives: 


A- - — A 

At ' 3 ~ dx 3 13 q 
< 94 ’ 

4 = ^ ? + V 

ft A.. 

A . . , - v -T q A . - V q A- 

ftAl 


(B.12) 


where Tk is a Cristoffel symbol defined in Eq.(B.17). With the above formulations, Eqs. 
(B.l)-(B.IO) can be written as follows 




(B.13) 
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(/*■<), + - rj^ + T’frpUiV" = -g 


where 


0 = U k k = ^ 

’ ax* 


• T ^ r, . /ac^ at/,- 2^ v 

+ s [(/• + «•) (a? + - 3 e »«) 

- <?’ K* + /-r)tf„g rj + 2rjg [(/* + w)£/.] 


-n ;.(.<■«) (g * g - - f »>,) 

) 


VTf^ + Atr) 


a?/, dUj 


+ ^lf- 2r ip^n - X05; 


ax? dxi 


+ 0 j> ( ^la _ r” r • - r* r- 1 

^ » l ax r **■ nj ,r ,n / 


+ rL^" 


and 


ff‘ J = 


ax* dx’ 

&T*a** ’ 


9ii = 


dX^_ dX ' k _ 
dx i dxi 


(B.14) 

(B.15) 

(B.16) 


here X k denotes the Cartesian coordinte system while x* represents a general coordinate 
system. The symbol TJ*, called the Christoffel symbol, is defined as 



dx % a 

dXPdXi 



(B.17) 


The equations for the turbulent kinetic energy, k = g'iuiUj = UiU and its dissipation rate e 
can be written as 


(P*0,t + 


dp£k 

ax* 


+ r in pkU n =gi 


a ' 

dx T 
+ Pk ~ pe 



dk 

dx n 


(B.18) 


(P £ ),t + 


dpUU 

dx* 


+ T\ nP eU n 



+ c, lfl ip t -c. 2 h ef +c , 3 «v|£|£ 

k k p dxj dx r 


de 

dx n 


(B.19) 
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where 


P k = -g J pu iUj \-^ + r kn U n 


2 , dUi dUj nT ^ nTT 2 

= ~z pkgi ‘ + » [sj? + a? - 2r ’> £/ ” - 3 effii 


(B-21) 


If we decompose P k into two parts, one due to the linear part of —pUiUj and the other due 
to the nonlinear part, then we may write 


p„ = pi 11 + pf 1 


(B.22) 


where 


P k = - o (p k + + g ki HT T-f + 


k iwr dU ^ dU i or n U 


p ( 2 ) _ /j fc J T . . f + r* U T 

Fk 9 13 l dx* + 1 mkU 


+ TUU" 


The nonlinear part of the cubic model, Eq. (B.7), in a general coordinate system is 


T . = Cl ^3^ + ^r*cP-r’i7 

•’ Cl C 2 V dx k dx> + 8x l -‘" U “ ■ 


r)TT k 


+ - fn lSy ) 

- [*“ (SS - - ww) - 

+ Cs l^ a? 9x> + <‘ u T * u 'dx> ■‘" p 3 n2ft v 

. r Ml [„« (WdUj_mr _ m9uz r , „ _ r , „ Widvz 

4 e 3 [ 9 Vdx*dx m dx‘ dx k dx l im 9 ikp dx m dx l 

„ „ at/ m at/,- ac/i ac/ m at;, ac/ m „ 

■ rP r? rr 77 4. _2 * _ •> p« 77 

ifc p 9 Qx l dx k dx m dx 1 dx k dx l tm q 

P mhdu - wr at* 

^■fc^P^m foJ +1 ifc 1 »mt> P t/g 9a .| + & r k0 x m i W t ' 

-0r?„t/,Ctr - ^nu,rsir + ryyW^Er 

I dUi rlJnrrr dUj q pmTrr 9t/, ^ pmr rr 

+i7»rj.r;t/,o' f tp) - |n, 9 «| 
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+ c/4 fSSS - d Sri,v P 4 + r 9U ‘ 


dx' dx l dxi dx' 


dxi 


qt dx 1 dxi 


, D dU l dU k dU k dU l dU k „ dU l 

9, w p dxi + foi dx 1 dx *' dxi k ‘ Up dx i 

, dU k dU l . „ at/* dU k dU k , 

+ qjU dx 1 dx' qi " 5x i + a** ax' 1 rjU 

flTjk flJT 

r’ u T‘ rj u t ir + - r* 


5x* 


^ x / * r 3 

dU k dU k , dU k „ , dU k k , 

SS 1 X T " V - -g^n/^U’ + ^r^tr 


- |n s9 «) 


where 


n > = •" (IP + (H - rLf/ ’) (f? + 


The scalar strain and rotation rates axe 


(B.25) 


(B.26) 

(B.27) 

(B.28) 


(s*) 2 = ^ 




ax* r ^ 9 


(S7*) 2 = 


■K£*w) 


f dU^_ 

dxi 


+r& 


y ') + (^ +r ^) © + r ^ 


'^ + r‘ 

dxi +i 'i 




dxi 

W" = Eq.{B. 39) 


W \ / du j * ,M 
fe7 + r ;^)(^r + r ^')J 


(B.29) 


(B.30) 

(B.31) 


B.3 Another form of the cubic model 


In terms of strain and rotation rates, the cubic Reynolds stress model can be written as 

pk 2 


2 pk 2 ( 1 \ 

pu{Uj = —p k gij + 0^ — 2 j -|- T+j 


(B.32) 
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where the nonlinear part, is 


Ty = 2A^ (fSi'S,, - i s< 2 ' S y) + 2A 2 ^- (fa, r si qi - in< 2 » 9ij ) 

+-A 3 ^2 9^ (^ipQqj ~ ^ipSqj^J “b 2^4/3— 9^9 ^Si p Sq T £l s j — QipS qT S s j^J 

-2A 5 ^- g pq g T3 Q, ip S qT £l 3j - -ClSQ, g {j + II S (s# ~ 9i?j 

-Is (B.33) 

where 



&Ui , _ 2r fc U 

dxi + dx«' ,J * 



\dx* 



© = S M , 5 (2) = gng^SprS'q , ft (2) = g^g^Sl^q 
OSH = g^g^g^SAq , J s = 0 , II S = \ (© 2 - S (2) ) 


and 


(B.34) 

(B.35) 

(B.36) 


(S ') 2 = !t*it's u s„ - ie 2 

(B.37) 

(fi') 2 = g^g'-a „,fi„ 

(B.38) 

O* C * Q* 

W* = g ij g kl g mn 

(S*) 3 

(B.39) 

where, 


Sij — Sij — -g 2 j 


Note that Eq. (B.33) appears to be more 
convenience for the CFD implementation. 

compact than Eq. (B.25) and may bring some 
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C Appendix: Equations in Cylindrical Coordinates 


C.l Mean equations 


Now, let us write all equations in a cylindrical coordinate system: x i = (x, r, 9). To accomplish 
this, we need to calculate the metric tensors g**, gij and the Christoffel symbol for 
cylindrical coordinates. Let X 1 = (x, y, z) be the cartesian system. The relation between the 
two systems is 

x = x , y = r cos9 , z = r sind (C.l) 


or 


x = x , 



6 — arctan (z/y) 


We may easily calculate 


dx l 

dXi 


' 1 

0 

0 ' 

dX { 

' 1 

0 

0 

0 

cos9 

sin9 


0 

cos9 

— r sin9 

.0 

—sin9/r 

cos9/r „ 

ox’ 

.0 

sin9 

r cos9 j 


(C.2) 


(C.3) 


The metric tensors g'i and g^ can then be obtained according to Eq. (B.16): 


r 1 

0 

0 ’ 


r 1 

0 

0 1 

0 

1 

0 

II 

0 

1 

0 

. 0 

0 

l/r 2 J 

.0 

0 

r 2 „ 


(C.4) 


and the Christoffel symbol r' fc can be obtained from Eq.(B.17) 


o 

o 

o 


o 

o 

o 


o 

o 

o 

0 0 0 

II 

f-4 

0 0 0 

II 

0 0 1/r 

o 

o 

o 

1 

o 

0 

1 


o 
1 — 4 

o 


(C.5) 


The contravariant velocity in the cylindrical coordinates is 


U' = ( U , V, Q) 


(C.6) 


where U and V are the axial and radial velocities, is the angular velocity. The corresponding 
covariant velocity can be obtained from 

Ui = 9ij U j = (U,V,r 2 n) (C.7) 
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With Eqs.(C.l)-(C.7), the equations for turbulent flows in a cylindrical coordinate system 
become 


Continuity equation 


dp dpU_ dpV dptt pV 
dt ^ dx ^ dr ^ dd ^ r 


= 0 


Momentum equations 


dpU dpU 2 dpUV dpUCl pUV 

I I I rvrt ^ 
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0 = "5 ^ “a ^ ~OQ ■* 

dx dr dd r 


(C.8) 


(C.9) 


(C.10) 


(C.ll) 


(C.12) 
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p=p + - z k 


k-E equations in Cylindrical coordinates 


dpk dpUk dpVk dpUk V _ d_ 

dt dx dr d9 r ^ dx 
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(C.13) 

(C.14) 

(C.15) 

(C.16) 

(C-17) 

(C-18) 
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C.2 Nonlinear part of turbulent stresses Tij 


After g ij , and for the cylindrical coordinate system are calculated, we may use Eq. 
(B.25) or Eq. (B.33) to calculate all the turbulent stresses automatically through a computer 
program. However, in the cylindrical coordinate system, most components of p t; , gij and 
are zero, therefore it is possiple to manually write down all the turbulent stresses to avoid 
many unnecessary null operations in the computer code. We write them here in a general 
form for the cubic model, so that model users can use their particular model coefficients for 
their applications. Note that with Eq. (A. 8), the coefficients C{ can be easily obtained from 
A { , or vice versa. 
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—2rfl^~ + V 2 + r 2 ft 2 - Vll 2 
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dd 

pk 4 \, 


1 

3 


+ C ^T 

£ 3 

+ c 4 

(C-19) 


T°° + T™ + ... + T% + Ttl - ^r 2 n : 

2 


7$ + 2$ + ... +2*? + 2*f - -r 2 n 3 


The terms . . . T${ in the above equations are listed below: 


The terms in t. 


T IXX 

19 


dU_ 

dx 


'dudu dudv_ + + 

dx dx + dr dx ^ dd dx ) 


dUdV d£<m 
dr y dx dr dr dr dd dr / 


rpxx 

20 


\_au_ (9U_dU_ dUdV 

+ ^~dd\dx dd + dr dd + 

= 0 



t 2 t = o 
T£= 0 

JTII J 1 IZ 


L 0*3 

2 * =0 
T|* = 0 

t 2 T = o 


'19 



t 2 t = o 
t% = o 

rpxx _ ft 
-^30 — u 

T iXX rpxx 

31 — 1 27 

Tg= 0 

T33* = 0 

T^=0 



dtfcw ac/ 

dx dx dr 
-dfi dU dr 2 QdV 


dV dU_dti\ dV ( dVdU 
dx ^ dd dx ) ^ dx \ dx dx 


dV dV dV_d& 
^ dr dx ^ dd dx 


dx dx 

dV dn 


T z -™ irx ~ tx + TV \ Tx ) 
!?«? = 0 
t 3 T=o 

T iXX rpxx 

40 “ -*36 

T ixx rpxx 

41 — ± Z1 

Tg = 0 


, dQ dfi \ 
dr dx dd dx ) 

{dti\ 2 


33 



34 


9U (dVdU_ 9V_dV_ dVdQ\ 
dr dr dr dr d9 dr J 


dv (du_du dUdv du_dn\ 

dr \dx dr + dr dr dd dr ) 


dv (dvw dvdv dvm\ 

dx y dx dr dr dr + d6 dr ) 



dx dr 

ov dci 



dCldU dr 2 Cl dV ~dCldCl\ 


^ h r ' 4 

dr dr d9 dr ) 

dvdci\ Tr dci 

H + rV — 

dr dx j dr 


rpxr 

-i 40 — 


+ 


dU (dUdU dU 
dr ydx dx dr 
dCl ( vdCldU dr 2 Cl dV 


dV dU_dCl\ dV_ ( 
dx dd dx J dr \ 


dVdU dVdV dV dti' 
dx dx ^ dr dx ^ dd dx , 


rpxr 

1 41 — 

rpxr 

1 42 ~ 


dr 

rpxx 

-*37 

«(r- 


+ 


dx dx dr dx 


+ r‘ 



,dCldU dr 2 Cl dV , dCldCl ' 


dx dx dr dx 


-f r 

dd dx 


T ta :t 
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rpxr 
-* 44 

rpxr 

rpxr 

1 46 
rpxr 

1 47 

T ixr 

48 

T ixr 
49 

T ixr 

50 

is 


dx dx 
ClfdUdU dVdV 


r \dx dd + 

Til 

0 

0 

0 

0 

0 

0 


dx dd 


+ r 



The terms in T„a : 



dfidC/ dr 2 CldV 
dx dx 
dCl dU 


+ r 




dr dx 
dr 2 Cl dV 


,dCldCl\ dU_{ 2 dCldU 
dd dx J + dr \ dx dr 


dr 2 Cl dV 2 

4 \- r 

dr dr 


+ r* 


,dCl dCl\ 


Tfo=rV 


Tg = 0 
T 2 t = 0 


d0 ' dr dd ' ' dd dd J 
dr r 2 dd dd) v 


'dU dCl dU_ 
dx dx + dr 


dU dV dUdV 1 8 XJ dV ' 
dx dx dr dr r 2 dd dd 


T. 1- 



dU dU 


dx dx 
dU dU 


du_dv_ du_du\ 

dr dx ^ dd dx ) 


drHl (dU_dU_ 
dr y dx dr 


dudv_ du_dci 

dr dr dd dr 


dU 


dx dd dr 


dV WdCl\ 
~dd + dddd) 


Tg = 0 

rpxB 

-* 25 


-rfi 


'dl/ dU dUdV dUdCl\ 
dx dr ^ dr dr ^ dd dr ) ^ 


VfdVdU dUdV_ dU_dCV 

7 Vd7 ~dd + ~d^~dd + ~dd~dd t 


dCl dCl 
~ddlfr 
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36 


(dU_dU_ dV dV 2 <mdtl\ 
\dx dr + dx dr + T dx dr J 


(dUdU dVdV dtldr 2 £l\ 
\dd!h + ltflb + ~dd~lfr- ) 


(dv du dvdv 

\ dx dx + dr dx + d9 dx ) 



The terms in t,,: 


ti; = 


+ 


+ 


dV (dU dV 
dx V dx dx 


1 dV ( dVdU 


dVdV dV_ an\ dV_ (dV^dU 
dr dx + dO dx ) + dr \ dx dr 


avdv dvm' 

^ dr dr ^ dO dr 


dO 


T T 2i f = -rft 

i s— “( 


\dx do 
dv an 


dVdV (TVdtt 
dr do + do dO 


y dx dx 
dV dU 


+ 


+ 

dV dn 

dr dr 


i_dvm" 
r 2 dO dO , 
dVdtl\ 


dVdV_ 

dx dO + dr dO + dO dO 


22 dO 

rprr rprr 

1 23 — 1 19 

rprr rprr 

^ 24 ~ a 20 

rrrrr rrrrr 

-*25 X 21 

rprr rprr 

-*26 ~ 1 22 



rprr 

1 29 — 


Q 2 V 

rrrrr T 

•*30 ~ „ 

T 

rprr rprr 

1 31 —-*27 

rprr rprr 

1 32 — 1 28 

rprr rprr 

x 33 — 1 29 

rprr rprr 

-*34 — -* 30 


~dV 

— n 2 — 

r 2 dO dr 

nvdv -dv 

1- n — 

n d0 dr 



dU dU_ dVdV 2 anan\ 
dr dx ^ dr dx ^ T dr dx ) 
dV 2 dQdQ\ 
dO +r dr dOJ 


dU_dU_ dV_ 
dr dO dr 


m n ~dv an _ _ . 

^ ~ 2rQ dr dr rV 

dCldU 

dx dr 




+ 


dV fdU dU dV dV dCt dr 2 Sl 


dr ^ dr dr dr dr dr dr 


+ 


+ r‘ 


an an\ 


- 


dr dO ) 



n (du du dv dv 


rriTT _ | v 

44 r \ dr dO 

rrtrr rprr 

^ 45 — 1 39 

__ ^ 2 dCl 

t:i - q 2 — 

46 dO 


,dCl dCl 


dr dO T dr dd 


) 


ft 2 . 


23 ? = 

/ 

rprr rprr 

J 48 — 1 44 

rprr rprr 

^ 49 “ 1 45 

rprr rprr 

a 50 ~~ 2 46 

T rrr rprr 

51 — -l 47 


■V 


47 


The terms in r.a : 


>dft dfi 


^ = dy / 2 a«dt/; a^nav; 2 

19 dx y dx dx dr dx dx 

_u dr 2 QdV 

+ r 2 dd \ dx d$ 

v /wan wan 

T *°- rV [tete + ~frfr + 



\ dVf 2 dU_dU 
)+ dr \ dx dr 


+ 


dr 2 Cl dV , dttdCl ' 


dr dr 


+ r* 


d6 dr 


t 2 i 


_ _o / 2 dn du 

r \ dx dO 


]_dV dfi\ _ Q [dVdV dVdV 1 dV < 
r 2 dO dO J y dx dx ^ dr dr ^ r 2 dO 


T£ = Cl 2 ^- - 

22 dO 

e _ 2 dn (dvdu 

23 dx \ dx dx 

dV dU_ 

r\ y-l i 


av% 

80 


+ 


dr 2 CldV 2 dCldn\ 

_) L p 2 I 

dr dO 80 80 j 


dV 



dr 

dV dV 
dx dO ' dr dO 


dV dydft\ 

dx "** 80 dx J 

dydfi\ 

dO 80 J 


drHlfdVdU dVdV_ 8V dCl' 
dr y dx dr dr dr dd dr 


rrrre „ / 2 dfi dn dr 2 ci dCl dd dfi\ 

■*24 - -rfi r ~a~~n~ + ~7T 5“ + ‘o7r"o7T 


dx dx dr dr dO dO ) 


^25 — ~rCl 


'dFd£7 dydy dFdf2\ y_(dV_dU dV dV dVdfA 
v dx dr + dr dr + d0 dr j + r i dx d0 + dr d0 + d0 d0 / 


■*26 — r i* dr S2 %0 

l dr d0 d0 dx 


FdVdfi 

+ r 90 90 


7^ = -rfi 3 - -y 2 
r 
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Til = 


rprO 

J 32 "" 

rprO 

J 33 “ 


n(dr 2 nav 2 anav\ vavan 


-n 


\ dr dd 
2 dr 2 Cl 


— r 


- nv an 

dd 

v 2 dv 


dd dr 

an 


+ 


r dd dd 


dr 

^ av v 2 av 2 av 
-nv— + ——-n 2 — 



dd 


dd 



du_au avav 2 anan\ av (au_du_ av_av_ 

dr dx dr dx dr dx ) ^ dd y dr dr dr dr 


+ 


an dr^cf 

dr dr 


dU dU av av 



dr dd dr dd 

avan avan\ andft 

dr dd dd dr ) r dr dd 


2 anau af^av 

dx dd + dr dd 
2 dV 

— n 2 — 


2 an an\ 

+ r dddd) 


r p rO 

^40 ~ 


d dd 

au_au_ avav 2 anan\ 

dd dx ^ dd dx ^ T dd dx J 


dV (dUdU dVdV dndr 2 n' 

+ ~dr + dd dr + dd dr 


+ 

t£ = 

rprO 

“M2 


an (dudu_ avav 

dr \dd~dd + dd dd 


2 an an\ 

+ r ~dd~dd) 


pxB 

■37 


v ( 2 anau ar 2 nav 2 

— r 1 1 - r 

r \ dx dr dr dr 

(dVdU d^&V dV 
\ dx dr ^ dr dr ^ dd 



t£ 

rprO 

■*44 

T£ 

rprO 
1 46 


rprv 

1 A 7 


Til 

rprO 
^ 49 

Til 

2 b ? 


av 


= -nv— + r 2 n 2 — + v 2 — 


an 



dr dr 

'au_au_ avav 

, 7777 + dr dd 


dr 


2 anan\ _ rQ (du_w 
+ r !h~aa) -r 


0 1 av_av an dr 2 n \ 

r \dr dr + dr dr + dr dr ) 


rpxO 

-M3 

, - a 

dd dr 

: rfi 3 + -V 2 

r 

n (duau avav 2 man\ 
: 7 \ dd dd + dd dd +r dd dd) 

rprO 
’ X 39 


rprO 

— 1 39 

Til 
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The terms in raa : 





, an at/ 9r 2 nav 


dx <9x 
dx 06 

on on 

Ox Ox 


+ 


dr Ox 


OrHlOV 
dr 06 
0r 2 il Oil 
Or Or 


2 0il0il\ Or 2 il 

+ T ae&) + 


Or 


t 


2 onou Or 2 nov 2 onon 

r 2 1 1 - r 

Ox Or Or Or 06 Or 


+ T d 



+ T‘ 




rpoe 

-*23 


= -rfi r - — — + 


Oil OU 0r 2 il OV 


Ox Or 


Or Or 

OV 

= -rW^ + r 2 ^— + V 

Or 



+ r 

, on 

06 


— — ) -ril (r 2 — 

06 06 ) V dx 

f\ + -{ 

Or J r \ 


OV OrHl OV 
Ox Or Or 


2 0il OV 

+ r 2 

06 06 


OilOU 0r 2 il OV , Oil Oil' 


06 


-HV 


Ox 06 


Or 06 


+ T 

06 06 


OV 

06 


+ ( V 2 + r 2 il 2 ) 


on 


rpOO 

1 29 



dr ' '86 

V 3 


OU OU OV OV 


06 Ox 
OU OU 


+ 


+ r 2 


onon\ 

06 Ox ' 06 Ox J 

OV OV Oil 0r 2 il 
06 Or + 06 Or + 06 Or 


/ 

\ Oil (dU_dU_ OV 
) + 06 \ 06 06 + ~06 


OV 

06 


+ r 2 


, Oil Oil' 
~06~06 


rpOO 

1 Z 7 — 


—2ril 


OV Oil 
~06~06 


+ rV 



rpOO 

-*38“ 


„ (OVOU OVOV 
' rfi ( Ox 06 + Or ~06 + 


OV 

~06 


rpOO 

-* 39 

rpOO __ 
“* 40 

K = 

l£ = 
T% = 


n ..dv ,„,an ,,,da 

-HV 1- r 2 il 2 h V 2 — 

06 06 06 

rptOO 

-*36 

r pOO 

-*38 

rpOO 

-*39 


™\ + V ( r l™ 

06 J r \ Ox 


OilOU Or 2 il OV 2 0il0il' 

06 + ~n r ~06 JrT ~06~06 
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^^vfeudu, 

44 r l 89 89 


* I V V V w 

-rO — — 


8U_8U_ 
89 dr 


8V8V 2 8080\ 
~89~8e +r 'de 89) 
8V8V 80, 8r 2 0 \ 
^ 89 dr + 89 dr J 


rn$$ rp$$ 

-^45 — -*39 

2- - _ 0 y + r 2 Q 2 ^ + V 2 — 

146 - iiV { dr + 89) + 8r +V 89 

T$ = — + 2 r0 2 V 

T 

rp$0 rrtOB 

^ 48 “ -*44 

rp99 rjiOO 

J 49 — ^39 

rp06 rpOS 

1 B0 — 1 A6 

r pG O rp99 

i 51 — X 47 


Other scalar quantities : 


n 1= 2^- 


80 0 8V V80\ _ _ 2 V 2 

8r + r 89 + r 00 ) + r 2 


(C.20) 


„ (8U\ (dvy 2 

n2= (a^J + (^J +r 

2 / at r \ 2 


'an' 

dx 


+ 


/acA 2 f&v 
\ dr J \ dr 


\ 2 80 8r 2 0 nr , 2 

) + a7^r +2ft 


+- 




w; + U*J +r \*>) - 2rn ^ + 2ry ^ + y2 


n 3 = n< 1} + n^ 2) + n 3 


( 3 ) 


(C.21) 

(C.22) 


nl 1} = 


'8U\ 3 8U (8U8V 8U_80\ dy_(dV_dU_ 8V 
k 8x J dx \ dr dx + 89 dx ) + dx \dx dx + dr 


8V8V_ 8V_80: 

8x 89 8x 


80 2 


'80 8U_ 808V 80 80\ „ (80\‘ 

8x 8x + dr 8x + 89 8x J + \ 8x ) 


ni 2) = 


8U ( 8U 8U 8U8V 8£80 
dr l 8x dr dr dr 89 dr 


yi( 


8V8U_ 8V 8V 8V 80' 
8x dr dr dr 89 dr 


(C.23) 


80 , 
+ aT 


'808U_ 808V_ 

dx dr dr dr 


+ + n( 


+rV 


'80' 

dr 


+ 2V0 


IllV = ^ 


_ 2 80 V -2 

+o 2 — + -o 2 

89 r 
8U 
89 


80 / 

a7 + rn ( 


80 80 
~89~dr ) 
80 8U 


8U8U_ dVdV j 
dr 89 ^ dr 89 ^ 


80 

dr 


80' 

89 


dodv 

8x dr ^ dr dr ^ 89 dr J 


80 80\ 


(C.24) 


fdUdU dU_8V 8U_dO\ 8V f 
\dx~89 + !fr~89 + ~89~d9 ) + 89 \ 


8V8U 8V8V 8V80' 
8x~89 + ~dr~89 + ~d0~d9 
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dn r2 /dndU MdV d£ldn\ (dUdU dVdV 
' dO T \dx dO + dr dO + dO dO ) ~ T \dOlfr + ~dOlfr + T ~dO~dr 


( dUdU avdv r2 Mdn\ (my _ av m 

r \ ae de do do do do) + \do) vu do +6V de 

+rV (OLHL + d J. d -L\ 

V dx dO dr dO ) l dx dO dr dO I 


-r 2 Vn^Y + r 2 ti 


,dV V 3 


(S') 2 = 7; 2 





2[H^J + {te) +r [to] + {d7j +2 [d^J 
r 2 f—Y 2—— 2 — — 2 —— 1(—Y 1 ( dV ' 

+ \ dr ) dr dx + dO dx dO dr + r 2 \ dO ) r 2 y dO 

+2 ^V +4 ^ +2 g.^ 


(*rj> = i f f » 2 + ,»« » + f »y + ^ («y + + 4si 2 _ 2 #^ 

2 \dx J dx dx y dr J y dr J dr dr dx 

n dUdU n dVdCl 1 (dU\ 2 1 (dV\ 2 SldV 1 

~ 2 dO dx 2 dO ~dr + r 2 [~dO ) + r 2 \dO ) ~ 4 7~d0 

W* = Eq.(CAl) (( 


C.3 Another form of Tij 


In terms of 5 tJ and f 2,j, the components of r tJ - can be written as 

T xx = 2 Ai^2~ (SllSll + 5 i 2 S 2 i + ^ 5i3«?3i — -S (2) ) 

+2A 2 ^- ^fi 12 fl 21 -f ^2^13^31 — + 2-43^- ^S 12 ft 21 + — 5x3^31^ 

+2-^4 25ii ^5 12 f2 2 x + — 5x3^31^ 4- 25 12 ^<S 22 f2 2 x 4 — -jS 2 3fi3x^ 

4" ^2 ‘-*13 ^3 2 fi 21 + ^533^31^1 

— 2-^-5 fll2 ^<S 22 fI 2 x 4- ^S^sfisx^ 4- ^3^13 ^<$32^21 4" — 

~ 3^^^ H s ^11 — 3®) — Is (s n S n 4- <Sx 2 5 2 i 4- — Sxs^si — — (C.29) 
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- — 


, P k 3 
+ A$ — y 
e 2 

_ pfc 4 

+2A4 — — 


2 -Ai^j- ^ 5 'ii 5 i 2 + 5 i 2522 + 5 x 3 5 3 2 ^ + 2 ^ 2 ^T^^ 13^32 

^12 ($11 — ^22) + "o’ (^ 13^32 ” ^IS^) 
r " 5 J 

5 n ( 5 u ni 2 + ^■ 5 'l 3 ^ 32 ^ + 5 x 2 ^ 523^32 
+ — 5 i 3 ( 5 3 lO U + ^ 5 33 fi 3 2 ^ ” ^12 ^‘S , 22 ‘S , 22 + ^ 52 3 5 3 2 ^ 

— ^^i 3 ^‘S' 3 i 5 i 2 + S32S22 + ^ 5 3 353 2 )j 
— 2 As^j- fii 2 (s 2 ifti 2 + 5 23 ^ 32 ^ + ^^13 ^ 31^12 + ^ 5 33 fi 32 ^ 

+/ 755 j 2 — Is (j>llSl 2 + S\ 3 .S 22 + ^ 5 i 3 5 3 2 ^ j 


= 2Ai^- (j>nSi3 + 5x2523 + ^^13 533^ + 2A2^~2~£}l2^: 
+^ 3 ^" 5 h «13 “ ^ 125 23 + ^ 2^23 ( 5 x 2 “ 5 3 2 )| 


+ 2A4 


pk 4 


5ix (5 n fii3 + 5x2^23 ) + 5x2(521^13 + $22^2z) 

~\~~^$ 13$ 32^23 — ^12 (52l5i3 + 522523^52 3 5 3 3^ 

^2 ^ 13 (^32 523 + ^ 5^33 533 ^| 


-2 A\ 


pk 4 


^12 (521^13 + 522^23) H 2 *^13 (531^13 + 532^23) 


+IIsS \2 — Is (5ii5i3 + 5x2523 + ^2 ^l 3 ^ 33 ^ 


= 2 a£( 


S 21 S 12 + 522522 H ^ 52 3 5 3 2 


- \ sm ) 


+ ^ 2^22 (531^x2 + ^-S , 33 fi 32 ^ 


3^32 — 3^ 2 ^) ^ ^2~5~ (^ 2 1^ 12 + ~^2$22' 

— 5i3f232 
T l 

j + 2S22 (^21^X2 + ~2 “-*23^32^ 

3^32^ 


— 5x3^32^ 

+ ^^23 (i>31^12 + -^33^32^ 

: e ) ( 

, ^ ^ 

5 2 i5i2 + S 22 S 22 “1 — j52 3 5 3 2 — - 
, r z 3 


‘32 


(C.30) 


(C.31) 


(C.32) 
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T T0 — ^521‘5i3 + S 22 S 23 A ^ 823833 ^ A 2A 2 ^-£l2ltt 

+-^3^“ (s 21 Sl 13 A S 22 Q 23 — Q 21 S 13 — ^^ 23833 ^ 


+2A 


pfc 4 


•S'?! (•S’llf^ 13 A <?12^23) + ‘S , 22(£21^13 + £ 22 ^ 23 ) 


+ ^2'S' 23 *S' 3 1^13 — ^21 (SnSu + «?12'^23^5'i3S , 33^ 


— ^2^23 ^£31 ‘ 


$13 4 7833833 


-2 A] 


pk 4 


)] 


^12 (Sufiis A S^f^) H 2^ 23 (^31^13 A 532^23) 


+//5S23 — I5 (5 2 i5i3 + S 22 S 23 A ^2'S , 23533^| 




— 2Aj^j- (s 31 S 13 A S 32 S 23 A ^'S'33‘S , 33 - — r 2 S^ 2 ^ 


+2A2~2~ ^ 31 J2i3 + ^32^23 ~ — r 2 fl^ 2 ^ + 2A3^y (531^13 + 532 ^ 23 ) 


+2A4-p- 2S > 3 x (*?ii^2i3 A S 12 Q 23 ) A 25 3 2 (SAfli 3 + S 22 Q 23 ) 


-2A, 


— ^2*^33 (<Sl3^31 + 523^32) 

pk 4 


^31 (*^11^13 A $12^23) A f232 (5 21 fii3 A 522^23) — -r 2 flSft 

O 

+ 1 Is ^‘S'33 — “T -2 ©^ — Is ($31813 A S 32 S 23 A -^2833833 — -r 2 S^ 2 ^ 


The scalars that appear in the above equations are as follows 


© — <?n + S22 H — 0833 
r 2 


5 (2) — 5ii5u + S 22822 H — 7833833 A 25i 2 5i2 H — -(5i 3 5i3 + S 22 S 23 ) 

T * V* 


2 

= 21712^12 H 2 (^13^31 + ^ 23 ^ 32 ) 


fiiSfi — fti 2 n 2 l(5n + ^ 22 ) + ^12^31 ^2 (*5*23 + *$32) + ^21^32“j(*S , X3 + S 3 1 ) 

+ -fi 3 1 fi 2 3(S 12 + S21) + ^fil3^3l(Sn + — S 33 ) + ^^23^32(^22 + ^^33) 


J s = 0, I/5 = ^(© 2 - S< 2 >) 


(C.33) 


(C-34) 

(C.35) 

(C.36) 

(C.37) 

(C.38) 

(C.39) 
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Two other scalars (S*) 2 and (SI*) 2 can be expressed as 


(s *) 2 = s (2) - ~e 2 , (si *) 2 = si (2) 

3 


(C.40) 


W* = 


(s*) z (^I 2 5 i 2 + + r 2 ^ 13^1 

+s ; 2 (sjjSi + 5 ^s 2 \ + i ) 
+^Si*3 (sj,Su + sja + is^sj.) 
+Sn (Si\S.‘2 + s; 2 s; 2 + is^) 
+s 2 - 2 (s 2 -,s; 2 + S 2 - 2 S 22 + ±s;,s ; 2 ) 
+isj 3 (sj,s; 2 + sj 2 s 2 - 2 + is; 3 s 3 - 2 ) 
(si'A + sia + is; 3 si) 
+is; 2 + 5 22 s 23 + 

+±s; 3 (r n s- a + s; 2 s 2 - 3 + is 3 ' 3 s; s ) 


(C.41) 


where 


S'u = 5u - ^0, s; 2 = S 22 - ^0, 5 3 *3 = S 33 - ir 2 ©, 


S* 12 = S 12 , 


$13 — $13i $23 ~ $ 23 


(C.42) 


Finally, the six components of Sij and the three components of Sl,j (note that 5,j = Sj x and 
Sljj = — SXjj) are 


Q _dU 1 (dU dV 

11 dx ’ 12 2 \ dr dx 


— o 5 ^ — <• 

ux 1 

a c _ 1 

522 - a7’ 323 

^ 1 (du dv 

SI 12 — r I ~x ^ 


„ 1 (dU 2 dn 

13 ~ 2 ( de + r ax 


„ 1 (dv 2 d(l\ n 2 dtl ir 

S23 = 2\-d9 +r fr)’ S33 ~ r ~de +rV 


n 1 feu 2 asi 

Ql3 - 2 ( dd r dx 


1 (dV 


^23 — T *rT 


2 v aa 


dr 2 Q, 

dr 
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D Appendix: Equations for Axisymmetric Flows 


Continuity equation 


dp 8pU dpV 
dt dx + dr 



Momentum equations 


P-1) 


dpU dpU 2 dpUV pUV 


dt 


dx 

d_ 
dx 
1 


dr r 

JdU 1 > 

2 (" + «•) &- 3®, 


dP_ 

dx 

d_ 

dr 


^ ( dU dV \ dr xx dr xr 

+ 7 ( " + ( a7 + to j + IT + IT + ; T " 


. /ay ay' 

+ «■> V a? + ^ 

1 


dpy dp£/F dpt/ 2 2 pV 2 

+ -4 — + -4 rptf + 


dt 


dx 


dr 


d_ 

dx 


. ( dv du' 

(p + Pt) f — f 


dP_ 

dr 


dx dr 


+ 2(p + Pt) 


dUL _ v 

Ur 3 , 

+ ; ( M + to) (f£ - |e) - £(». + w) (V - ir 2 e) 


+ 


dr,, dr,, 


+ 


dx dr 


o T 00 + 


dpr 2 n dpr 2 rtu dpr 2 nv d 

dt dx dr dx 


d_ 

dr 


(p + Pr) 


dr 2 f2 

dr 


(p + Pr) 


dr 2 tt 

dx 


1, % dr 2 ffc 

+ -(a+w)- gr 


- - jr [(/» + ar)r 2 n] + ^ 
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